Run Modular Discretization

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Introduction

The objective of this notebook is to discretize HRUs in a watershed and create a weighting matrix that describes the hillslope-stream network connectivity in a catchment. This notebook is meant to replace the spatial discretization workflow (DynatopSpatialFunctionExplicitReaches) for the modified dynamic TOPMODEL described by Mahoney et al. (2022) J. Hydrol.

To run this notebook, it is necessary to first install kernels for Python and R and install arcpy to the environment from which the notebook is run (see instructions here). Instances when switching the kernel is necessary will be denoted in a Markdown cell.

Step 1: Load in necessary libraries for both R and Python

We will be using both R and python in this example. Python will be called via the reticulate package. We will be using the python installation from arcgis pro. A major advantage of using the R Notebook document over, for example, Jupyter Notebook is the ability to store variables generated by both python and R in the environment. Thus we can use variables generated in R with the arcgis API, have the user-friendliness of Rstudio, and have the replicability of a Notebook document.

setwd('C:/Users/david/OneDrive/Desktop/Mahoney Research/2-Mahoney-EPA/modeling-streamflow-permanence/')

# Load in libraries
library(reticulate)
library(sp)
library(raster)
library(dynatopmodel)
library(topmodel)
library(ecbtools)
library(cluster)
library(factoextra)
library(tidyverse)
library(NbClust)
library(basicClEval)
library(ggplot2)
library(mapview)
library(doParallel)
library(foreach)
use_python('C:/Program Files/ArcGIS/Pro/bin/Python/envs/arcgispro-py3/python.exe',required=T)

Warning in Sys.setlocale("LC_CTYPE", ctype) :
  using locale code page other than 65001 ("UTF-8") may cause problems

arcpy <- import("arcpy")
arcgis <- import("arcgis")
py_math <- import("math")
py_numpy <- import("numpy")
ipython <- import("IPython")

Step 2: Import Rasters and calculate TWI

Note: This code will calculate the TWI using the D8 flow direction. It is possible to import a DINF or MFD raster here as well and use this to run the analysis, assuming that it is already in the GDB.

# Read in the rasters that we'd need for HRU creation
dem <- raster('SpatialInputData/fr1meterDEM.tif')
soils <- raster('SpatialInputData/FRSoils1mMask.tif')
TWI.dinf <- raster('C:/Users/david/OneDrive/Desktop/Mahoney Research/2-Mahoney-EPA/modeling-streamflow-permanence/SpatialInputData/fr_twi_dinf.tif')
TWI.dinf[TWI.dinf==0] <- NA # We do this because the areas outside of the raster have a value 0 and we want these to be set to NA

Step 3: Create a HRU raster based on the TWI, soils, and other data using R

Now, we will switch over to R and create an HRU map of for the watershed. Make sure to switch the kernel from python (esri) to R

Step 3.1: Install libraries and read in raster files

Relevant libraries include spatial and raster libraries, cluster analyses, and tidyverse. Note, we are creating HRUs with the following rasters:

• DEM

• Soils

• TWI

Other relevant layers that could be used include:

• LULC

• Ksat

• Porosity

• Aspect

• Curvature

• Temperature

• Precipitation

Note: there is some degree of customizability here. Any relevant layer could be added.

Step 3.3: Create an elbow plot showing the number of HRUs which should be created

The total within sum of squares for each cluster is used to estimate the variability for each iteration.

Note: this code takes a decent amount of time to run.

wss.raster <- function(k) {
    # Run the K-means analysis
    kmeans.raster <- raster.kmeans(raster.stack, k, iter.max=10, nstart = 10, geo = T, geo.weight = 1)
    
    # Load in the data and get into vector form
    HRU.vector <- data.frame('hru' = as.vector(kmeans.raster))
    TWI.vector <- data.frame('TWI' = as.vector(TWI.dinf.cut))
    soils.vector <- data.frame('soils' = as.vector(soils))
    elev.bands.vector <- data.frame('Elev.bands' = as.vector(elev.bands))
    
    # Put the raster vectors into a cleaned data frame
    rasters.df <- data.frame(HRU.vector,TWI.vector,soils.vector,elev.bands.vector)
    rasters.df.clean <- rasters.df[complete.cases(rasters.df),]
  
    # Calculate the within cluster sum of squares
    within.ss <- wcss(rasters.df.clean[,-1],rasters.df.clean[,1])
    
    # Return the total within cluster sum of squares
    return(within.ss$WCSS)
}

set.seed(123)

# Compute and plot within sum squares for k = 1 to k = 40
k.values <- 1:30

# extract within sum of squares for 1-40 clusters
wss_values <- map_dbl(k.values,wss.raster)

Warning: Quick-TRANSfer stage steps exceeded maximum (= 20786250)
Warning: Quick-TRANSfer stage steps exceeded maximum (= 20786250)

# plot within sum of squares vs cluster (HRU) number
plot(k.values,wss_values,type='b',pch=19,frame=F,xlab='Number of clusters',ylab='Total within-cluster sum of sq')

Step 3.4: Plot the clusters for reference

We will look at the various clusters created and how they appear within relation to one another

# Load in the data and get into vector form
HRU.vector <- data.frame('hru' = as.vector(HRU))
TWI.vector <- data.frame('TWI' = as.vector(TWI.dinf.cut))
TWI.nc.vector <- data.frame('TWI.nc'=as.vector(TWI.dinf))
soils.vector <- data.frame('soils' = as.vector(soils))
elev.bands.vector <- data.frame('Elev.bands' = as.vector(elev.bands))
elev.bands.nc.vector <- data.frame('Elev.bands.nc' = as.vector(dem))

# Put the raster vectors into a cleaned data frame
rasters.df <- data.frame(HRU.vector,TWI.vector,TWI.nc.vector,soils.vector,elev.bands.vector,elev.bands.nc.vector)
rasters.df.clean <- rasters.df[complete.cases(rasters.df),]

# Plot the elev bands vs twi with HRU as color
ggplot(data=rasters.df.clean, aes(x=TWI.nc,y=Elev.bands.nc,col=hru))+geom_point()

# Note - the straight lines are due to the cuts that we creted. We could just as easily not include these cuts and other geometries would come through

Step 3.5: Sort the HRUs by TWI and add in the stream raster

In this step, we will first sort HRUs by TWI value (high TWI = high number HRU), then we will “burn in” a stream network that has been previously delineated. This represents the geomorphic stream channel, to which the watershed uplands will contribute water.

The raster need to have the following characteristics:

• The resolution, extent, and snapping should be the same as the DEM

• The raster is derived from a Flow accumulation raster with an upstream contributing area threshold applied

• The threshold should create a stream network that is similar to the stream mapped in the field

• The reaches should be created using the stream link/order tools in ArcGIS

• The shapefile created from these stream tools should be buffered such that a contiguous stream network raster can be created

• Here, a buffer of 2 cell widths was used (~1.524 * 2 m)

## Calculate the mean TWI for each HRU - relabel the HRUs such that high valued TWI is the highest number HRU and low valued TWI is the lowest HRU
## Note - this isn't necessary for dynatopmodel, but may be necessary for dynatop
mean.twi <- data.frame(matrix(nrow=length(unique(HRU)),ncol=5))
names(mean.twi) <- c('HRU.name','atb.bar.unsort','id.new','atb.bar','id')
mean.twi$HRU.name <- unique(HRU)

hillslope.hru.matrix <- as.matrix(HRU)
TWI.dinf.matrix <- as.matrix(TWI.dinf)

mean.twi$atb.bar.unsort <- zonal(x=TWI.dinf,z=HRU,fun='mean')[,2]
mean.twi$atb.bar <- sort(mean.twi$atb.bar.unsort,decreasing=F)
mean.twi$id.new <- match(mean.twi$atb.bar.unsort,mean.twi$atb.bar)
mean.twi$id <- 1:length(unique(HRU))

reclassify.matrix <- data.frame(matrix(nrow=length(unique(HRU)),ncol=2))
names(reclassify.matrix) <- c('is','becomes')
reclassify.matrix$is <- mean.twi$HRU.name
reclassify.matrix$becomes <- mean.twi$id.new

HRU <- reclassify(HRU,reclassify.matrix)

## Overlay the stream network atop the Hru raster
stream <- raster('SpatialInputData/FR_stream_network.tif')
stream[stream==0] <- NA
HRU <- HRU*10+max(max(unique(stream)),100)               # Uniquely identifying the upland HRUs
HRU[stream>0] <- stream[stream>0]                     # Identify the cells from the HRU raster that overlap with the stream, set those equal to the reach number 
writeRaster(HRU,"SpatialInputData/FR_HRU.tif",overwrite=T)

image(HRU)

twi.zonal <- zonal(x=TWI.dinf,z=HRU,fun='mean')
TWI.zonal <- reclassify(HRU,twi.zonal)

image(TWI.zonal)

stream.lump <- stream
stream.lump[stream >= 1] <- 1
HRU.lump <- HRU
HRU.lump[stream.lump>0] <- stream.lump[stream.lump>0]

writeRaster(HRU.lump,"SpatialInputData/FR_HRU_lump.tif",overwrite=T)
writeRaster(stream.lump,"SpatialInputData/FR_stream_lump.tif",overwrite=T)

image(HRU.lump)

Step 3.6 Visualize the HRU/Stream Raster

Note - I’m going to do this with the mapview function in R, but this could be done using the arcpy and arcgis modules in python. The reason for not using python is because if we interrupt the kernel we’ll have to reload variables from above

mapView(HRU, maxpixels =  2102000) + mapView(TWI.zonal, maxpixels =  2102000) + 
  mapView(HRU.lump,maxpixels =  2102000) + mapView(stream.lump, maxpixels =  2102000) +
  mapView(stream,maxpixels =  2102000)

Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition

Step 4: Create a flow weighting raster

In this step, we will create the weighting matrix to run the dynatopmodel analysis. We do this twice. Once for the “Explicit matrix” which shows the connectivity of upland cells that contribute to a specific reach. A second time for the “Lumped Matrix” which shows the connectivity of cells to a “lumped” stream network. As of 8/2/2022, the “dynatopmodel” package is not compatible with explicit reaches.

Step 4.1: Calculate the explicit downstream matrix

This process returns the value of the HRU immediately downstream of a cell. This algorithm strictly works for a D8 flow direction, but could be expanded for DInf at a later time.

Make sure to switch over to the R kernel

# First, create the explicit matrix
library(sp)
library(raster)

HRU <- raster('SpatialInputData/FR_HRU.tif')
fdr <- raster('SpatialInputData/fr_fdr')


HRU.matrix <- as.matrix(HRU)
fdr.matrix <- as.matrix(fdr)

fdr.downstream.matrix <- matrix(nrow=nrow(HRU.matrix),ncol=ncol(HRU.matrix))

for (i in 1:nrow(HRU.matrix)) {
  for (j in 1:ncol(HRU.matrix)) {
    if (is.na(fdr.matrix[i,j])) {
      fdr.downstream.matrix[i,j] <- NA
    } else {
      if (fdr.matrix[i,j]==1) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i,j+1]
      } else if (fdr.matrix[i,j]==2) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i+1,j+1]
      } else if (fdr.matrix[i,j]==4) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i+1,j]
      } else if (fdr.matrix[i,j]==8) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i+1,j-1]
      } else if (fdr.matrix[i,j]==16) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i,j-1]
      } else if (fdr.matrix[i,j]==32) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i-1,j-1]
      } else if (fdr.matrix[i,j]==64) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i-1,j]
      } else if (fdr.matrix[i,j]==128) {
        fdr.downstream.matrix[i,j] <- HRU.matrix[i-1,j+1]
      } else {
        fdr.downstream.matrix[i,j] <- NA
      }
    }
  }
}

Step 4.2: Calculate the explicit weighting matrix which shows the number of cells from one HRU to contribute to downstream HRUs.

In this step, we will create a weighting matrix that outputs the proportion of flow generated in one cell that will be redistributed to downstream cells.

This works by:

• identifying all cells that belong to a certain HRU

• identifying all cells which belong to a downstream HRU

• identifying the cells where both of these conditions are true

• returning the number of cells where both are true

In this regard, the output will be all cells which belong to HRU i that flow into HRU j.

Note: I rewrote some of this code. The run time goes from about 50 minutes to approximately 50 seconds.

print('Running weighting matrix...')

[1] "Running weighting matrix..."

no.HRU.explicit <- length(unique(HRU))
HRU.explicit.list <- unique(HRU)

downstream.weighting.matrix <- matrix(nrow=no.HRU.explicit,ncol=no.HRU.explicit)
start.time <- Sys.time()

for (i in 1:no.HRU.explicit) {
  print(HRU.explicit.list[i])
  for (j in 1:no.HRU.explicit) {
    downstream.weighting.matrix[i,j] <- sum((HRU.matrix==HRU.explicit.list[i])&(fdr.downstream.matrix==HRU.explicit.list[j]),na.rm=T)
  }
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
[1] 7
[1] 8
[1] 9
[1] 10
[1] 11
[1] 12
[1] 13
[1] 14
[1] 15
[1] 16
[1] 17
[1] 18
[1] 19
[1] 20
[1] 21
[1] 22
[1] 23
[1] 24
[1] 25
[1] 26
[1] 27
[1] 28
[1] 29
[1] 30
[1] 31
[1] 32
[1] 33
[1] 110
[1] 120
[1] 130
[1] 140
[1] 150
[1] 160
[1] 170
[1] 180
[1] 190
[1] 200
[1] 210
[1] 220
[1] 230
[1] 240
[1] 250
[1] 260
[1] 270
[1] 280
[1] 290
[1] 300

end.time <- Sys.time()
end.time-start.time

Time difference of 53.847 secs

options("scipen"=100, "digits"=5)

rownames(downstream.weighting.matrix) <- unique(HRU)
colnames(downstream.weighting.matrix) <- unique(HRU)

weighting.matrix <- downstream.weighting.matrix/rowSums(downstream.weighting.matrix)

rownames(weighting.matrix) <- unique(HRU)
colnames(weighting.matrix) <- unique(HRU)

rowSums(weighting.matrix)

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
 21  22  23  24  25  26  27  28  29  30  31  32  33 110 120 130 140 150 160 170 
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
180 190 200 210 220 230 240 250 260 270 280 290 300 
  1   1   1   1   1   1   1   1   1   1   1   1   1 

Step 4.3: Initialize the lumped downstream matrix

Note that the current implementation of dynatopmodel does not support the explicit streams, as far as I am aware. Dynatop may, but I haven’t investigated the model in depth as of 8/3/2022

# Next, create the lumped matrix

HRU.lump <- raster('SpatialInputData/FR_HRU_lump.tif')
fdr <- raster('SpatialInputData/fr_fdr')

HRU.lump.matrix <- as.matrix(HRU.lump)
fdr.matrix <- as.matrix(fdr)

fdr.downstream.lump.matrix <- matrix(nrow=nrow(HRU.lump.matrix),ncol=ncol(HRU.lump.matrix))

for (i in 1:nrow(HRU.lump.matrix)) {
  for (j in 1:ncol(HRU.lump.matrix)) {
    if (is.na(fdr.matrix[i,j])) {
      fdr.downstream.lump.matrix[i,j] <- NA
    } else {
      if (fdr.matrix[i,j]==1) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i,j+1]
      } else if (fdr.matrix[i,j]==2) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i+1,j+1]
      } else if (fdr.matrix[i,j]==4) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i+1,j]
      } else if (fdr.matrix[i,j]==8) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i+1,j-1]
      } else if (fdr.matrix[i,j]==16) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i,j-1]
      } else if (fdr.matrix[i,j]==32) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i-1,j-1]
      } else if (fdr.matrix[i,j]==64) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i-1,j]
      } else if (fdr.matrix[i,j]==128) {
        fdr.downstream.lump.matrix[i,j] <- HRU.lump.matrix[i-1,j+1]
      } else {
        fdr.downstream.lump.matrix[i,j] <- NA
      }
    }
  }
}

Step 4.4: Calculate the lumped weighting matrix which shows the number of cells from one HRU to contribute to downstream HRUs.

Note: I have rewritten the code here. The previous run time for this chunk was 10 minutes. Now it is about 10 seconds.

print('Running weighting matrix...')

[1] "Running weighting matrix..."

no.HRU.lump <- length(unique(HRU.lump))
list.HRU.lump <- unique(HRU.lump)

downstream.weighting.matrix.lump <- matrix(nrow=length(unique(HRU.lump)),ncol=length(unique(HRU.lump)))
start.time <- Sys.time()
for (i in 1:no.HRU.lump) {
  print(list.HRU.lump[i])
  for (j in 1:no.HRU.lump) {
    downstream.weighting.matrix.lump[i,j] <- sum((HRU.lump.matrix==list.HRU.lump[i])&(fdr.downstream.lump.matrix==list.HRU.lump[j]),
                                                 na.rm=T)
  }
}

[1] 1
[1] 110
[1] 120
[1] 130
[1] 140
[1] 150
[1] 160
[1] 170
[1] 180
[1] 190
[1] 200
[1] 210
[1] 220
[1] 230
[1] 240
[1] 250
[1] 260
[1] 270
[1] 280
[1] 290
[1] 300

end.time <- Sys.time()
end.time-start.time

Time difference of 10.993 secs

options("scipen"=100, "digits"=5)

rownames(downstream.weighting.matrix.lump) <- unique(HRU.lump)
colnames(downstream.weighting.matrix.lump) <- unique(HRU.lump)

weighting.matrix.lump <- downstream.weighting.matrix.lump/rowSums(downstream.weighting.matrix.lump)

rownames(weighting.matrix.lump) <- unique(HRU.lump)
colnames(weighting.matrix.lump) <- unique(HRU.lump)

rowSums(weighting.matrix.lump)

  1 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 
  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
300 
  1 

Step 5: Create the Groups Matrix

Next, we will create a “Groups” matrix which will specify the initial parameters for dynamic TOPMODEL.

Note - I did notice some differences between atb.bar from the 20 HRUs discretized here and the 2- HRUs discretized with dynatopmodel. It is related to how the HRUs are discretized. There is more variability with the dynatopmodel approach.

Items we will need to include in the table:

• id - the ID of the channel, typically channels are ranked 0-100 and HRUs are ranked based on TWI with ID > 100

• tag - same as id

• chan.no - the channel number (note - only for IDs classified as channels)

• area_pc - the percent area that the HRU takes up

• area - the area m^2 that the HRU takes up

• sbar - the average slope of the HRU (m/m)

• atb.bar - the average TWI for the HRU

• gauge.id - the rain gauge contributing to the HRU

• catch.id - the catchment in which the HRU is situated

• srz_max - maximum depth of root zone parameter value initialized as 0.1

• ln_t0 - ln t0 parameter value initialized as 7

• m - exponential decline of transmissivity parameter initialized as 0.01

• srz0 - initial root zone storage parameter initialized as 0

• td - td parameter initialized as 1

• vchan - v chan parameter (for channels only) parameterized as 1000

• vof - overland v parameter (for HRUs only) parameterized as 100

• k0 - k0 parameter initialized as 1e+8

• CD - CD parameter (unused) initialized as 0.1

• sd_max - sd_max parameter initialized as 0.5

• pe_fact - potential evapotranspiration factor used to scale PET, initialized as 1

• vof_fact - vof factor initialized as 1

• rain_fact - rain factor used to scale precip amount, initialized as 1

• mann.n - mannings n initialized as 0.01

• S0 - initial slope, parameterized as 0.1

• ex_max - ex max parameterized as 1

Step 5.1 Create the groups matrix for the explicit reaches

# Enter the cell size 
cell.size <- as.numeric(readline('Enter the resolution of the raster in m (for 10-m raster enter 10): '))
1.524

# Create functions to determine the area of the raster
area_cell <- function(raster.matrix,raster_ID) {
  area_c <- sum(raster.matrix==raster_ID, na.rm=T)
  return(area_c)
}

# Create function to determine the zonal mean (unused)
mean_raster <- function(raster.matrix, HRU.matrix, raster_ID) {
  mean_rast <- mean(raster.matrix[HRU.matrix==raster_ID],na.rm=T)
  return(mean_rast)
}

# Calculate the slope 
slope <- terrain(dem,opt='slope',units='tangent')

# Initial ize the groups 
groups.explicit <- data.frame(matrix(nrow=length(unique(HRU)),ncol=26))
names(groups.explicit) <- c('id','tag','chan.no','order','area_pc','area','sbar','atb.bar','gauge.id','catch.id','srz_max','ln_t0','m','srz0','td','vchan','vof','k0','CD','sd_max','pe_fact','vof_fact','rain_fact','mann.n','S0','ex_max')

# Assign the groups
groups.explicit$id <- unique(HRU)
groups.explicit$tag <- unique(HRU)
groups.explicit$chan.no[1:length(unique(stream))] <- unique(stream)
groups.explicit$order <- 1:length(unique(HRU))
groups.explicit$area_pc <- sapply(X=groups.explicit$id,FUN=area_cell,raster.matrix=HRU.matrix)/sum(!is.na(HRU.matrix))*100
groups.explicit$area <- groups.explicit$area_pc*sum(!is.na(HRU.matrix))*cell.size^2/100
groups.explicit$sbar <- zonal(x=slope,z=HRU,fun='mean')[,2]
groups.explicit$atb.bar[(1+length(unique(stream))):length(unique(HRU))] <- zonal(x=TWI.dinf,z=HRU,fun='mean')[(1+length(unique(stream))):length(unique(HRU)),2]
groups.explicit$gauge.id <- 1
groups.explicit$catch.id <- 1
groups.explicit$srz_max <- 0.1
groups.explicit$ln_t0 <- 7
groups.explicit$m <- 0.01
groups.explicit$srz0 <- 0
groups.explicit$td <- 1
groups.explicit$vchan[1:length(unique(stream))] <- 1000
groups.explicit$vof[(1+length(unique(stream))):length(unique(HRU))] <- 100
groups.explicit$k0 <- 1e+8
groups.explicit$CD <- 0.1
groups.explicit$sd_max <- 0.5
groups.explicit$pe_fact <- 1
groups.explicit$vof_fact <- 1
groups.explicit$rain_fact <- 1
groups.explicit$mann.n <- 0.01
groups.explicit$S0 <- 0.1
groups.explicit$ex_max <- 1

Step 5.2: Create the groups matrix for the lumped reaches

# Calculate the slope 
slope <- terrain(dem,opt='slope',units='tangent')

# Initialize the groups 
groups.lump <- data.frame(matrix(nrow=no.HRU.lump,ncol=26))
names(groups.lump) <- c('id','tag','chan.no','order','area_pc','area','sbar','atb.bar','gauge.id','catch.id','srz_max','ln_t0','m','srz0','td','vchan','vof','k0','CD','sd_max','pe_fact','vof_fact','rain_fact','mann.n','S0','ex_max')

# Assign the groups
groups.lump$id <- list.HRU.lump
groups.lump$tag <- list.HRU.lump
groups.lump$chan.no[1:length(unique(stream.lump))] <- unique(stream.lump)
groups.lump$order <- 1:no.HRU.lump
groups.lump$area_pc <- sapply(X=groups.lump$id,
                         FUN=area_cell,raster.matrix=HRU.lump.matrix)/sum(!is.na(HRU.lump.matrix))*100
groups.lump$area <- groups.lump$area_pc*sum(!is.na(HRU.lump.matrix))*cell.size^2/100
groups.lump$sbar <- zonal(x=slope,z=HRU.lump,fun='mean')[,2]
groups.lump$atb.bar[(1+length(unique(stream.lump))):length(unique(HRU.lump))] <- zonal(x=TWI.dinf,z=HRU.lump,fun='mean')[(1+length(unique(stream.lump))):length(unique(HRU.lump)),2]
groups.lump$gauge.id <- 1
groups.lump$catch.id <- 1
groups.lump$srz_max <- 0.1
groups.lump$ln_t0 <- 7
groups.lump$m <- 0.01
groups.lump$srz0 <- 0
groups.lump$td <- 1
groups.lump$vchan[1:length(unique(stream.lump))] <- 1000
groups.lump$vof[(1+length(unique(stream.lump))):length(unique(HRU.lump))] <- 100
groups.lump$k0 <- 1e+8
groups.lump$CD <- 0.1
groups.lump$sd_max <- 0.5
groups.lump$pe_fact <- 1
groups.lump$vof_fact <- 1
groups.lump$rain_fact <- 1
groups.lump$mann.n <- 0.01
groups.lump$S0 <- 0.1
groups.lump$ex_max <- 1

Step 6: Routing Matrix

In this step, we will create routing matrix that bins lists the length of each cell to the outlet and the frequency.

# Determine the flow length from each cell in the stream to the catchment outlet
outFlowLength <- arcpy$sa$FlowLength('SpatialInputData/fr_fdr',"DOWNSTREAM")
# Note - the next line only needs to be saved once. 
#outFlowLength$save(paste0(getwd(),'/SpatialInputData/Flow_length.tif'))

Flow.length <- raster('SpatialInputData/Flow_length.tif')
Flow.length.stream <- stream.lump
Flow.length.stream[!is.na(stream.lump)] <- Flow.length[stream.lump==1]
writeRaster(Flow.length.stream,'SpatialInputData/Flow_length_stream.tif', overwrite =T)

input.units <- readline('Was the original DEM in units of ft or m (enter ft or m): ')
ft

breaks.hist <- 10      # User defined number of breaks

RoutingTable <- data.frame(matrix(nrow=breaks.hist,ncol=2))
names(RoutingTable) <- c('flow.len', 'prop')

if(input.units == 'ft') {
  Flow.length.stream.fix <- Flow.length.stream * 0.3048
  Flow.length.vector <- na.omit(as.numeric(as.vector(Flow.length.stream.fix)))
  
  RoutingHist <- hist(Flow.length.vector,breaks=seq(min(Flow.length.vector),max(Flow.length.vector),length.out=breaks.hist+1))
  RoutingTable$flow.len <- RoutingHist$breaks[-1]
  RoutingTable$prop <- RoutingHist$counts/sum(RoutingHist$counts) 
  
} else if (input.units == 'm') {
  Flow.length.vector <- na.omit(as.numeric(as.vector(Flow.length.stream)))
  RoutingHist <- hist(Flow.length.vector,breaks=seq(min(Flow.length.vector),max(Flow.length.vector),length.out=breaks.hist+1))
  RoutingTable$flow.len <- RoutingHist$breaks[-1]
  RoutingTable$prop <- RoutingHist$counts/sum(RoutingHist$counts)
} else {
  'invalid response, RoutingTable not calculated.'
}

Step 7: Compile spatial data and compare original Spatial Function to outputs from this script.

We will create a list that compThe original DynatopSpatialFunctionExplicitReaches function outputs the following information:

• DEM

  • Contains DEM raster name

• Soils

  • Contains Soils raster name

• layers - contains all layers to make HRUs

  • DEM

  • Soils

  • TWI (named atb)

• disc - contains all elements to make the weighting matrix and groups table from the lumped parameterization

  • groups - the groups matrix

  • layers - same as above

  • chans - multiband raster with 1) channel labels (just equal to 1 for lumped) and 2) channel props, which describes the proportions of the cell which is occupied by the channel (here this can just be set equal to one if a 1.524-m raster is being used).

  • cuts - the number of cuts made to the layers

  • area.thresh - the area threshold below which HRUs will be lumped with other HRUs. This isn’t utilized here and can be set equal to 0.

  • hru - the HRU raster generated from the discretize process

  • weights - the weighting matrix

• RoutingTable - a routing table showing the flow distance to the outlet and the proportion of channel cells that fall within the distance bin - this had 30 cuts

• explicit.disc - contains all elements to make the weighting matrix and groups table from the explicit parameterization

  • groups - explicit groups matrix

  • layers - should be more or less the same as disc

  • chans - This time should be explicit with respect to Reach ID

  • cuts - should be the same as disc

  • area.thresh - should be same as disc

  • hru - This will be different than disc bc there should be explicit representations for each reach

  • weights - explicit weighting matrix

• explicit.RoutingTable - same as above, this had 5 cuts, however

• explicit.ChanTable - this is a table containing all of the reach data, including: “link_no,” “ds_link,” “us_link_1,” “us_link_2,” “stream_order,” “stream_length_ft,” “stream_drop_ft,” “stream_slope_ftpft,” “stream_straightline_length_ft,” “US_area_m2,” “width_m2”

# DEM
DEM <- dem

# Soils
Soils <- soils

# Layers 
layers <- raster.stack

# disc
disc <- list()
disc$groups <- groups.lump
disc$layers <- layers
disc$chans <- addLayer(stream.lump,stream.lump)
names(disc$chans) <- c('chans','chanprops')
disc$cuts <- data.frame(matrix(nrow=1,ncol=length(names(layers))))
colnames(disc$cuts) <- names(layers)
for (i in names(disc$cuts)) {print(i);disc$cuts[1,i] <- length(unique(layers[[i]]))}

[1] "TWI.dinf.cut"
[1] "soils"
[1] "elev.bands"

disc$cuts <- as.matrix(disc$cuts)
disc$area.thresh <- 0
disc$hru <- HRU.lump
weighting.matrix.lump.in <- weighting.matrix.lump
colnames(weighting.matrix.lump.in)[1] <- 'R'
rownames(weighting.matrix.lump.in)[1] <- 'R'
disc$weights <- as.matrix(weighting.matrix.lump.in)

# RoutingTable
RoutingTable <- RoutingTable

# explicit.disc
explicit.disc <- list()
explicit.disc$groups <- groups.explicit
explicit.disc$layers <- layers
explicit.disc$chans <- addLayer(stream,stream.lump)
names(explicit.disc$chans) <- c('chans','chanprops')
explicit.disc$cuts <- disc$cuts
explicit.disc$area.thresh <- 0
explicit.disc$hru <- HRU
weighting.matrix.explicit.in <- weighting.matrix
names.explicit.matrix <- paste0('R',unique(stream))
rownames(weighting.matrix.explicit.in)[1:length(unique(stream))] <- names.explicit.matrix
colnames(weighting.matrix.explicit.in)[1:length(unique(stream))] <- names.explicit.matrix
explicit.disc$weights <- as.matrix(weighting.matrix.explicit.in)

# explicit.RoutingTable (unused)
explicit.RoutingTable <- RoutingTable

# explicit.ChanTable
drn <- shapefile('SpatialInputData/FR4000StreamNet.shp')
bkf.width <- 12.16*(drn$USAreaM2/1000^2*0.386102)^0.42*0.3048  # estimate the width of the channel using the regional curves from Ashland Berry's thesis. Note this isn't perfect bc this equation is being extrapolated to smaller channels.
explicit.ChanTable <- as.data.frame(cbind('link_no'= drn$LINKNO,'ds_link'= drn$DSLINKNO,
                                    'us_link_1' = drn$USLINKNO1,'us_link_2'=drn$USLINKNO2,
                                    'stream_order'= drn$strmOrder,
                                    'stream_length_m' =drn$length_m,
                                    'stream_length_ft'=drn$Length,
                                    'stream_drop_ft'=drn$strmDrop,
                                    'stream_slope_ftpft'=drn$Slope,
                                    'stream_straightline_length_ft'=drn$StraightL,
                                    'US_area_m2'=drn$USAreaM2,
                                    'width_m2' = bkf.width))

# Compile the final spatial table
dynatop.spatial <- list('DEM'=DEM,'Soils'=Soils,'layers'=layers,'disc'=disc,
                        'RoutingTable'=RoutingTable,'explicit.disc'=explicit.disc,
                        'explicit.RoutingTable'=explicit.RoutingTable,
                        'explicit.ChanTable'=explicit.ChanTable)

Appendix A: Unused code

This chunk is a scratch workspace for old code

tot.cells <- data.frame(matrix(nrow=no.HRU.lump,ncol=2))
names(tot.cells) <- c('tot.cells','tot.cells.2')

for (j in 1:no.HRU.lump) {
  print(list.HRU[j])
  tot.cells$tot.cells[j] <- sum((HRU.lump.matrix==list.HRU[1])&(fdr.downstream.lump.matrix==list.HRU[j]),na.rm=T)
  #tot.cells$tot.cells.2[j] <- length(which((HRU.lump.matrix==unique(HRU.lump)[1])&(fdr.downstream.matrix==unique(HRU.lump)[j])))
  
}

parallel::detectCores()
n.cores <- parallel::detectCores() - 1
#create the cluster
my.cluster <- parallel::makeCluster(
  n.cores, 
  type = "PSOCK"
  )
#check cluster definition (optional)
print(my.cluster)

#register it to be used by %dopar%
doParallel::registerDoParallel(cl = my.cluster)

#check if it is registered (optional)
foreach::getDoParRegistered()

get.ds.cells <- function(i,HRU.matrix,HRU,fdr.downstream.matrix) {
  col <- matrix(nrow=1,length(unique(HRU)))
  for (j in 1:length(unique(HRU))) {
    temp.col <- length(which((HRU.matrix==unique(HRU)[i])&(fdr.downstream.matrix==unique(HRU)[j])))
    col[j] <- temp.col
  }
  return(col)
}

start.time <- Sys.time()
downstream.weighting.matrix.test <- foreach(i=1:length(unique(HRU)),
                                       .combine =rbind) %dopar% {
                                         temp.col <- get.ds.cells(i,HRU.matrix,HRU,fdr.downstream.matrix)
                                       }

end.time <- Sys.time()
end.time-start.time

stopCluster(cl = my.cluster)

options("scipen"=100, "digits"=5)

rownames(downstream.weighting.matrix.test) <- unique(HRU)
colnames(downstream.weighting.matrix.test) <- unique(HRU)

weighting.matrix.test <- downstream.weighting.matrix.test/rowSums(downstream.weighting.matrix.test)

rownames(weighting.matrix.test) <- unique(HRU)
colnames(weighting.matrix.test) <- unique(HRU)